function C = jacobian_f_gps(x0)

% your code here
eps = 1e-3;
%C = zeros(length(x0));

for j=1:numel(x0)
   x = x0;
   x(j)=x0(j)+eps;
   y2=f_gps(x, 0, 0);
   x(j)=x0(j)-eps;
   y1=f_gps(x, 0, 0);
   
   C(:,j)=(y2-y1)/(2*eps);
end
